Noise-induced phase transitions in neuronal networks 
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Using an exactly solvable cortical model of a neuronal network, we show that, by increasing the 
intensity of shot noise (flow of random spikes bombarding neurons), the network undergoes first- and 
second-order non-equilibrium phase transitions. We study the nature of the transitions, bursts and 
avalanches of neuronal activity. Saddle-node and supercritical Hopf bifurcations are the mechanisms 
of emergence of sustained network oscillations. We show that the network stimulated by shot noise 
behaves similar to the Morris-Lecar model of a biological neuron stimulated by an applied current. 

PACS numbers: 87.19.lj, 87.19.1n, 87.19.1c, 05.70.Fh 



In the brain, interactions among neurons lead to di- 
verse collective phenomena such as phase transitions, 
self-organization, avalanches, and brain rhythms |l|, |2|. 
A non-equilibrium second-order phase transition was ob- 
served in human bimanual coordination 3-5]. Stimula- 
tion of living neural networks by electric fields causes a 
first-order phase transition [6| . There are evidences that 
brain rhythms, epileptic seizures, and the ultraslow os- 
cillations of BOLD fMRI patterns also emerge as a result 
of non-equilibrium phase transitions pj. The Hopf and 
saddle-node bifurcations are generic mechanisms for the 
emergence of oscillations in nonlinear differential equa- 
tion models [8]. These mechanisms were found in the 
Morris-Lecar model of a biological neuron [9|, [lOj . In the 
context of brain rhythms, the Hopf bifurcations were dis- 
cussed within mean-field cortical models [7] and for net- 
works of randomly connected integrate-and-fire neurons 
Ill4l4| . At the present time, understanding of nature of 
collective phenomena in the brain is elusive. For a com- 
plete description of a phase transition it is not enough to 
identify the bifurcation. It is also necessary to find criti- 
cal phenomena that accompany it. In statistical physics, 
exactly solved models largely help us to understand phase 
transitions and critical phenomena [15| . 

In the present paper, we propose an exactly solvable 
cortical model with stochastic excitatory and inhibitory 
neurons stimulated by shot noise (a flow of random spikes 
bombarding neurons). We show that shot noise stimu- 
lates first- and second-order non-equilibrium phase tran- 
sitions in collective dynamics of neuronal networks. The 
first-order phase transition occurs as a discontinuous 
transition from low to high neuronal activity. Avalanches 
precede the transition. The saddle-node and supercritical 
Hopf bifurcations are the mechanism of emergence of sus- 
tained network oscillations. We find the order parameter 
for the continuous phase transitions and critical phenom- 
ena that accompany them in activity fluctuations. We 
show that the model exhibits collective excitability simi- 
lar to excitability of the Morris-Lecar neuron stimulated 
by an applied current. 

Model. We study a cortical model composed of N e ex- 
citatory and N{ inhibitory neurons. N e + Ni = N is the 



network size, g e u) = N e /j\/N is the fraction of excitatory 
(inhibitory) neurons. Neurons are randomly connected 
with probability c/N and form a classical random graph 
with Poisson degree distribution and the mean degree 
c. The network is locally tree-like and has the small- 
world properties |17H19J | similar to ones found in brain 
networks [20(. To take into account noise playing an im- 
portant role in the brain dynamics 2l423J|. we assume 



that neurons are bombarded by random spikes repre- 
sented by Dirac delta functions, I{t) — J2iQ^(^ ~ U)i 
where ti are arrival times of spikes and q is their am- 
plitude. This kind of random input is a so-called shot 
noise. Random spikes represent spontaneous release of 
neurotransmitters in synapses (synaptic noise) and ran- 
dom spikes arriving from the remote part of the cortex. 
According to Schottky's result, in the case of the Pois- 
son distribution of interspike intervals, the power spectral 
density S(lu) is proportional to the mean frequency ui sn 
of spikes, S(ui) — 2q 2 u> sn . We assume that the probabil- 
ity to receive £ random spikes during the integration time 
t is Gaussian, G(£) = G exp[— (£ — (n)) 2 /2a 2 ] where a 2 
is the variance, (n) — w sn T is the mean number of spikes 
arriving during the time interval r, and Go is the normal- 
ization constant, J^TLo @(0 = 1- We will use (n) as the 
control parameter characterizing the shot noise intensity. 
Neurons also receive delta-like spikes from active neigh- 
bors. The spikes mediate interaction among neurons. 
The total input I(t) includes spikes from shot noise and 
excitatory and inhibitory presynaptic neurons. We define 
the input V n to a neuron with index n, n = 1, 2, . . .N, as 
integral of I{t) over the time interval [t — r, t] . This gives 



V n (t)=Zq + kJe + Ui.. 



(1) 



where £, k, and I are the numbers of spikes arriving dur- 
ing the time interval [t—T, t] from shot noise, active presy- 
naptic excitatory and inhibitory neurons, respectively. 
We assume that efficacies of synaptic connections with 
excitatory and inhibitory neurons are uniform and equal 
to J e (J e > 0) and J t (Jj < 0), respectively The num- 
bers k and I are random and are determined by activity 
of presynaptic neurons during the interval [t — r, t] ■ The 
network structure is encoded in the adjacency matrix. 



In our model, neurons are tonic and the firing fre- 
quency f(V) versus input V is the Heaviside function 
f{V) = fQ(V - V th ) where V th is a threshold. / is 
the same for both excitatory and inhibitory neurons. If 
Jt < 1 and spike emission times of neurons are uncorre- 
cted, then during the time interval [t — r, t] , each active 
presynaptic neuron contributes to V n (t) either one spike 
with probability rf or none with probability 1 — rf. 

We consider stochastic neurons like those of 2J-|27( . It 
means that response of neurons to an input is a stochas- 
tic process. Such a stochastic behavior might be caused 
by cellular noise and intensive bombardment by random 
spikes. Two rules determines dynamics of the cortical 
model: (1) If the input V n (t) at an inactive excitatory 
(inhibitory) neuron n at time t is at least a certain thresh- 
old Vth, then this neuron is activated with probability 
fi e r (/J-iT) and fires spikes. (2) An active excitatory (in- 
hibitory) neuron n is inactivated with probability p e T 
(pir) if V n (t) < Vth- In our model, l//i e and l/pi are 
of the order of the first-spike latencies of excitatory and 
inhibitory neurons (from 6 to 100 ms, in the cortex). 
We introduce a parameter a = pi/p e - If a < 1, in- 
hibitory neurons have a larger response time than exci- 
tatory neurons, a may be both larger and smaller than 
1. We also introduce a dimensionless activation threshold 
il = Vth/Je- ^ is of the order of 15-30 in living neuronal 
networks and about 30 — 400 in the brain. 

Rate equations. Behavior of the model is described 
by the fractions p e (t) and pi(t) of active excitatory and 
inhibitory neurons, respectively, at time t. We will call 
them 'activities'. We assume that activities are changed 
slightly during the integration time r. Using the rules 
formulated above and the methods developed in [26l - l28j . 
in the limit N — >■ oo, we find explicit rate equations, 

pe{t) =Fe(t)[l-p e (t)]-MePe(*) + Me*e(Pe(t),/Oi(*)), 

pi{t) = Fi(t)[l-pi(t)} - mpiit) + Hi<iri(p e {t),pi(t)), (2) 

where p = dp/dt. ^ e U){Pe, Pi) is the probability that, at 
given activities p e and pi, the input to a randomly chosen 
excitatory (inhibitory) neuron is at least Q. F e and F( 
represent fields acting on excitatory and inhibitory neu- 
rons. In the case of the classical random graph, we find 

®i(Pe,Pi) = ^e(Pe,Pi) = ^{Pe,Pi) = 
EZ,^0 Q ( kJ e + l^Hq-^Je)G(OPk(gePe^Pl(g t Pr^,^) 

where c = erf, Pk(g e PeC) and Pi(giPiC) are the probabil- 
ities that, during the time interval r, a randomly chosen 
neuron receives k spikes from excitatory and I spikes from 
inhibitory neurons, respectively. Here Pk(c) = c k e~ c /kl. 
Algorithm. In simulations, we builded a directed net- 
work, linking neurons with the probability c/N. We di- 
vided time into intervals of width At = r. At each time 
step, for each neuron we calculated input Eq. (J]) given 
that each active presynaptic neuron contributes a spike 
with probability rf. The number of random spikes (shot 




FIG. 1. Points 1,2, and 3 represent solutions of the steady 
state equation p — ^(p, p) for the cases (n) < n c \ (solid line), 
n c i < (n) < n C 2 (dashed line), and (n) > n C 2 (dotted line). 



noise) in this input was generated by the Gaussian pro- 
cess G(£). Then, we updated states of neurons, using the 
rules formulated above. In the paper, numerical calcu- 



lations are presented for parameters N — 10 



10 J , 



n = 30, rf = 0.1, / = p e , and g t = 0.25. We use 
1/ p e e 1 as time unit and J e = 1 as input unit. Fol- 
lowing [16j, we choose J b — — 3J e . We use q — J e and 
a 2 = 10 for the amplitude and variance of shot noise. 

Steady states. If F e = Fi = 0, then in a steady state, 
we have p e = pi = p and p is a solution of the steady 
state equation, p ~ *f>(p,p). A graphical solution of the 
equation is displayed in Fig. [T] If shot noise intensity (n) 
is either sufficiently small or sufficiently large, then there 
is only one solution, either point 1 or point 3. These fixed 
points correspond to the steady states with low and high 
neuronal activity, respectively. In an intermediate range 
n c \ < (n) < n C 2 there are three fixed points (1,2, and 3). 
At the critical point (n) = n c i, points 2 and 3 coalesce. 
Points 1 and 2 coalesce at (n) — n C 2- From Fig. [T] one 
sees that the coalescence occurs when d^{p,p)/dp = 1. 
Together with the steady state equation, this condition 
determines n c \ and n C 2- While the fixed points depend 
on (n), but not on a, their local stability depends on 
both (n) and a. It is determined by eigenvalues of the 
Jacobian of Eqs. ([2|), 



m 



-l + d^/dpe d^/dpi \ 
ad^/dpe -a + ad^/d Pi ) 



(4) 



calculated at the fixed points. The eigenvalues are 



A± = -\{Jn 



1 



J 22 )±|a/(Jh- J 22 ) 2 + 4J 12 J 21 , (5) 

where Jij are the entries of the Jacobian. If X± < at a 
fixed point, then this point is stable (attractor). If A± > 
0, then the point is unstable. If one of the eigenvalues 
A± is positive and the other is negative, then the point 
is saddle. If ReA± < and ImA± ^ 0, the point is a 
stable spiral. If ReA± > and ImA± ^ 0, the fixed point 
is an unstable spiral. The real and imaginary parts of 
A± determines the relaxation rate, "f r — — ReA + , and the 
frequency, 7$ = ImA + , of damped oscillations. 

Phase diagram. Analyzing the local stability of the 
fixed points in the a — (n) plane (see Table |T|) , we found 



TABLE I. Stability of the fixed points f , 2, and 3 in the 
regions f, If, and III on the phase diagram in Fig. [5] where 
s=stable, sd=saddle, u=unstable, sp=spiral, lc=limit cycle. 
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FIG. 2. (Color online) (n) — a plane of the phase diagram 
of the cortical model, (n) is the shot noise intensity, a is 
the ratio of the response time of excitatory neurons to the 
response time of inhibitory neurons. The phase regions, the 
phase boundaries, and the parameters used in numerical cal- 
culations are explained in the text, at ~ 0.80. Lines 1 and 2 
represent two scenarios discussed in the text. 



the phase diagram displayed in Fig. rjj In regions Ia- 
Ie, the network relaxes exponentially to the stable fixed 
point 1 (if perturbations are small) . In regions lb and Ha, 
relaxation to the stable fixed point 3 is exponential while, 
in regions Ic and lib, the relaxation occurs in the form of 
damped oscillations. In regions Ilia and Illb, the fixed 
point 3 is an unstable point surrounded by a limit cycle. 
These are the regions with sustained network oscillations 
about the point 3. In Fig. [2 the phase boundaries rep- 
resented by the dashed and solid lines are determined by 
the conditions Ji(p^) = and 7r(p < - 3 - > ) = 0, respectively. 
Vertical lines represent lines (n) = n c \ and n c2 . 

First- order phase transition. The pattern of collective 
behavior depends on a and (n). First we consider the 
case a > at (at is the critical value below which sus- 
tained oscillations appear) . In simulations and numerical 
solution of Eqs. @, we increased the noise level (n) from 
zero (region la) to a value in region Ha (or lib) above the 
critical point n c2 and afterwards decreased it again to a 
value below n C 2 (see line 1 in Fig. [5]). Neuronal activity 
undergoes a jump at (n) — n C 2 (n C 2 ~ 18.8 in Fig. (3Jb)). 

Avalanches. Near n c2 , we observed bursts of neuronal 
activity caused by avalanches (see the inset in Fig.(3][a)). 
We studied avalanches, analyzing spike time series by 
use of the standard method (see [291 ] or the recent work 
[30|]). The avalanche size distribution P(s) is represented 
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FIG. 3. (Color online) (a) Avalanche size distribution P(s) 
versus size s at (n) = 18.8. Inset: temporal activity of excita- 
tory neurons near the first-order phase transition. Time t is 
in unit l/fi e - (b) Hysteresis in neuronal activity for increasing 
and decreasing noise level (n). In simulations, a = 0.85. 



in Fig. [3ja). When (n) is close to n c2 , P(s) is powerlaw, 
P(s) oc s~ CT , in a broad range of s. Using the maximum 
likelihood estimate 3l|, we found a sw 1.53 and the cor- 
responding p-value is p = 0.89 (the closeness of p to 1 
shows that the fit is good). This estimation agrees with 
the experimental data [29j, |30( and the standard mean- 
field ex ponent a = 3/2 obtained for other exactly solved 



models [32H34J . Thus, in the cortical model, bursts and 



avalanches are precursors of the first-order phase tran- 
sition. Another mechanism of avalanches based on self- 
organized criticality is discussed in [2] . 

Hysteresis. If (n) decreases from a value above n C 2 to 
a value below n C 2 , the network activity remaines as high 
as it was above n c2 (see Fig. Gl^b)). The activity falls to 
a low value only at (n) = n c \ < n c2 . Hysteresis occurs 
for a and (n) in regions lb and Ic and is absent in regions 
Id and Ie where the fixed point 3 is unstable. Hysteresis 
was observed, for example, in living neural networks 35] 
and in simulations of thalamocortical systems 36] . 

An analytical analysis of the cortical model gives us 
deeper understanding of the first-order phase transition. 
Solving the equation p = ty(p,p) at (n) near n c2 and 
substituting p' 1 ) into Eq. ([5]), we find that, in regions lb 
and Ic the relaxation rate to the low activity state p^ is 
7 r oc (n C 2 — (n)) 1 ' 2 , i.e., 7 r — > when (n) —> n C 2- This 
phenomenon is known as critical slowing down. Solv- 
ing Eqs. <j2j) with weak white-noise forces F a (i) (F a (t) oc 
1/y/N), which mimic forces caused by finite-size effects, 
we found, using the linear perturbation analysis |37j, that 
the power spectral density (PSD) of the activity fluctu- 
ations, S(u>) — (dp a (u)Sp a (—u>)}, has a sharp zero fre- 
quency peak, S(ui) oc l/(u> + %)■ If ( n ) ^ n c2, the peak 
maximum increases as S ma x oc 1/7^ (note that in the 
high activity state, there is no singularity around n C 2). 
If size N is finite, then j r remains nonzero, though very 
small, even at the critical point due to finite-size effects 
that smooth the transition observed in simulations. This 
leads to a finite value of S ma x at ri c2 . 

N on- equilibrium, phase transitions. An interesting sce- 
nario of collective dynamics takes place if inhibitory neu- 
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FIG. 4. Network oscillations near (a) the saddle-node 

((n) = 18.805, n C 2 = 18.8) and (b) supercritical Hopf ((n) = 
34, n C 3 = 36) bifurcations, (c) Amplitude (solid line) and fre- 
quency (dashed line) of network oscillations versus (n) (from 
a numerical solution of Eqs. (H}). At (n) > n C 3, the oscil- 
lations are damped. Time t is in units 1//U. e , and a = 0.75. 



rons respond slower than excitatory neurons and a < a t . 
In our analysis, we used simulations of the model, ana- 
lytical and numerical solution of Eqs. ©. At a given a, 
we start from (n) — and increase the shot noise inten- 
sity (n) (see line 2 in Fig. [5]). The neuronal network goes 
from region la with the single fixed point 1 into region 
Id or Ie where its dynamics is determined by three fixed 
points (see Table [I]). At (n) = n C 2 the points 1 and 2 coa- 
lesce and the network undergoes a non-equilibrium phase 
transition due to the saddle- node bifurcation. Above n C 2, 
the system is in region Ilia (or Illb) with sustained os- 
cillations (an unstable spiral or unstable fixed point sur- 
rounded by a limit cycle). Then, at (n) — n C 3(a) on the 
boundary between regions Ilia and lib, the network un- 
dergoes a phase transition due to the supercritical Hopf 
bifurcation. Above n c ^ the network enters region lib 
with damped network oscillations (a stable spiral in the 
point 3). The bifurcations were observed in the Morris- 
Lecar model stimulated by an applied current when the 
I— V relation is N-shaped [9j. Thus, collective behav- 
ior of neuronal networks stimulated by a flow of random 
spikes bombarding neurons is similar to behavior of a 
single neuron stimulated by an applied current. From 
this similarity it follows that the network acts as an 'in- 
tegrator', when it is close to the saddle-node bifurcation, 
and as a 'resonator', when it is close to the Hopf bifurca- 
tion [lOj. The important difference between the models 
is that, in our case, we deal with collective dynamics and 
phase transitions and it is necessary to find the order 
parameter and describe critical phenomena. 

Saddle-node bifurcation. Above the saddle-node bifur- 
cation, (n) > n c2 , the network oscillations emerge with 
a large amplitude (see Fig.rjja)) and their frequency in- 
creases from zero aswa ((n) — n^) 1 ^ 2 (see Fig, life)). 
This behavior is generic for the bifurcation |8H10l|. We 
suggest that the frequency is the order parameter of the 
transition. In simulations, at (n) below n C 2, we observed 
random almost identical bursts of activity that are the 
manifestation of critical fluctuations (inset of Fig. [5£a)). 
The bursts differ from ones found near the first-order 
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FIG. 5. (Color online) (a) Variance of the reciprocal of inter- 
burst interval A versus (n) about the saddle-node bifurcation. 
Inset: bursts of activity and a single burst at (n) — 18.76. 
Time t is in units l//i e . (b) The peak maximum of the PSD 
S(u)) of fluctuations versus (n) above the supercritical Hopf 
bifurcation (n C 3 ~ 80.5). Inset: Temporal neuronal activity 
and S{uj) versus ui at (n) = 82.5. In simulations, a — 0.55. 



phase transition in Fig. [3£a). We assume that the bursts 
are non linear excitations generated by finite-size effects 
|38j . The variance of the reciprocal of the interburst in- 
terval A, ([A" 1 - (A" 1 )] 2 ), where (A" 1 ) is the mean 
burst frequency, has a peak at n c i (see Fig. Ufa)) and 



behaves as \(n) —n C 2\~ about n C 2- Burst properties and 
generation mechanisms need further investigations. 

Supercritical Hopf bifurcation. Near the supercriti- 
cal Hopf bifurcation at n C 3, network oscillations differ 
from oscillations near the saddle-node bifurcation (see 
Figs. 0Ja) andSJb)). A decrease of the oscillation am- 
plitude, A ex (n C 3 — (n)) 1 ' 2 , and the relaxation rate, 
7 r oc |n C 3 — (n)\, signals the supercritical Hopf bifur- 
cation (see Fig. Etc)). This behavior is generic for the 
bifurcation [8]. The amplitude is the order parameter 
for the transition. We found these results, expanding 
i&(p e ,Pi) in Eqs. @ in a series in Sp a (t) = p a (t) — p^ 
around the fixed point 3 and holding terms up to 0{5p i a ) 
inclusively. Then, we solved Eqs. $2$ in region Ilia, us- 
ing the averaging theory [8J . If (n) tends to n C 3 from the 
'paramagnetic' region lib with damped oscillations, then 
fluctuations of neuronal activity are increased, signaling 
the bifurcation. The PSD of the fluctuations has a res- 
onance peak at the frequency ujq of damped oscillations: 
S{uo)/S max w 4C 2 /[(1 ~ x 2 ) 2 + 4C 2 z 2 ], where x = w/uo, 
C = "f r /ujQ is the damping ratio, and the peak maximum 
Smax oc I/7 2 (see also [12rll4J ]). These results agree well 
with our simulations in Fig. Etb). They may help to 
understand the damping of alpha rhythms by visual or 
auditory stimuli [39|, |4fJ . 

In conclusion, using the exactly solvable cortical model 
of neuronal networks with stochastic neurons, we showed 
that a flow of random spikes bombarding neurons con- 
trols behavior of the network. By increasing the flow, 
network oscillations appear at the saddle-node bifurca- 
tions and disappear at the supercritical Hopf bifurcation. 
Burst, avalanches, and critical fluctuations of neuronal 
activity are precursors of the non-equilibrium transitions. 
Within the model, neuronal networks are excitable sys- 



terns having dynamical behavior similar to the Morris- 
Lecar model of a biological neuron. 
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